{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as npl\n",
    "import scipy.linalg as spl\n",
    "import sklearn.metrics as sklm\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [6, 6]\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [],
   "source": [
    "irisD = {'virginica': np.array([[6.3, 3.3, 6. , 2.5],\n",
    "       [5.8, 2.7, 5.1, 1.9],\n",
    "       [7.1, 3. , 5.9, 2.1],\n",
    "       [6.3, 2.9, 5.6, 1.8],\n",
    "       [6.5, 3. , 5.8, 2.2],\n",
    "       [7.6, 3. , 6.6, 2.1],\n",
    "       [4.9, 2.5, 4.5, 1.7],\n",
    "       [7.3, 2.9, 6.3, 1.8],\n",
    "       [6.7, 2.5, 5.8, 1.8],\n",
    "       [7.2, 3.6, 6.1, 2.5],\n",
    "       [6.5, 3.2, 5.1, 2. ],\n",
    "       [6.4, 2.7, 5.3, 1.9],\n",
    "       [6.8, 3. , 5.5, 2.1],\n",
    "       [5.7, 2.5, 5. , 2. ],\n",
    "       [5.8, 2.8, 5.1, 2.4],\n",
    "       [6.4, 3.2, 5.3, 2.3],\n",
    "       [6.5, 3. , 5.5, 1.8],\n",
    "       [7.7, 3.8, 6.7, 2.2],\n",
    "       [7.7, 2.6, 6.9, 2.3],\n",
    "       [6. , 2.2, 5. , 1.5],\n",
    "       [6.9, 3.2, 5.7, 2.3],\n",
    "       [5.6, 2.8, 4.9, 2. ],\n",
    "       [7.7, 2.8, 6.7, 2. ],\n",
    "       [6.3, 2.7, 4.9, 1.8],\n",
    "       [6.7, 3.3, 5.7, 2.1],\n",
    "       [7.2, 3.2, 6. , 1.8],\n",
    "       [6.2, 2.8, 4.8, 1.8],\n",
    "       [6.1, 3. , 4.9, 1.8],\n",
    "       [6.4, 2.8, 5.6, 2.1],\n",
    "       [7.2, 3. , 5.8, 1.6],\n",
    "       [7.4, 2.8, 6.1, 1.9],\n",
    "       [7.9, 3.8, 6.4, 2. ],\n",
    "       [6.4, 2.8, 5.6, 2.2],\n",
    "       [6.3, 2.8, 5.1, 1.5],\n",
    "       [6.1, 2.6, 5.6, 1.4],\n",
    "       [7.7, 3. , 6.1, 2.3],\n",
    "       [6.3, 3.4, 5.6, 2.4],\n",
    "       [6.4, 3.1, 5.5, 1.8],\n",
    "       [6. , 3. , 4.8, 1.8],\n",
    "       [6.9, 3.1, 5.4, 2.1],\n",
    "       [6.7, 3.1, 5.6, 2.4],\n",
    "       [6.9, 3.1, 5.1, 2.3],\n",
    "       [5.8, 2.7, 5.1, 1.9],\n",
    "       [6.8, 3.2, 5.9, 2.3],\n",
    "       [6.7, 3.3, 5.7, 2.5],\n",
    "       [6.7, 3. , 5.2, 2.3],\n",
    "       [6.3, 2.5, 5. , 1.9],\n",
    "       [6.5, 3. , 5.2, 2. ],\n",
    "       [6.2, 3.4, 5.4, 2.3],\n",
    "       [5.9, 3. , 5.1, 1.8]]), 'setosa': np.array([[5.1, 3.5, 1.4, 0.2],\n",
    "       [4.9, 3. , 1.4, 0.2],\n",
    "       [4.7, 3.2, 1.3, 0.2],\n",
    "       [4.6, 3.1, 1.5, 0.2],\n",
    "       [5. , 3.6, 1.4, 0.2],\n",
    "       [5.4, 3.9, 1.7, 0.4],\n",
    "       [4.6, 3.4, 1.4, 0.3],\n",
    "       [5. , 3.4, 1.5, 0.2],\n",
    "       [4.4, 2.9, 1.4, 0.2],\n",
    "       [4.9, 3.1, 1.5, 0.1],\n",
    "       [5.4, 3.7, 1.5, 0.2],\n",
    "       [4.8, 3.4, 1.6, 0.2],\n",
    "       [4.8, 3. , 1.4, 0.1],\n",
    "       [4.3, 3. , 1.1, 0.1],\n",
    "       [5.8, 4. , 1.2, 0.2],\n",
    "       [5.7, 4.4, 1.5, 0.4],\n",
    "       [5.4, 3.9, 1.3, 0.4],\n",
    "       [5.1, 3.5, 1.4, 0.3],\n",
    "       [5.7, 3.8, 1.7, 0.3],\n",
    "       [5.1, 3.8, 1.5, 0.3],\n",
    "       [5.4, 3.4, 1.7, 0.2],\n",
    "       [5.1, 3.7, 1.5, 0.4],\n",
    "       [4.6, 3.6, 1. , 0.2],\n",
    "       [5.1, 3.3, 1.7, 0.5],\n",
    "       [4.8, 3.4, 1.9, 0.2],\n",
    "       [5. , 3. , 1.6, 0.2],\n",
    "       [5. , 3.4, 1.6, 0.4],\n",
    "       [5.2, 3.5, 1.5, 0.2],\n",
    "       [5.2, 3.4, 1.4, 0.2],\n",
    "       [4.7, 3.2, 1.6, 0.2],\n",
    "       [4.8, 3.1, 1.6, 0.2],\n",
    "       [5.4, 3.4, 1.5, 0.4],\n",
    "       [5.2, 4.1, 1.5, 0.1],\n",
    "       [5.5, 4.2, 1.4, 0.2],\n",
    "       [4.9, 3.1, 1.5, 0.2],\n",
    "       [5. , 3.2, 1.2, 0.2],\n",
    "       [5.5, 3.5, 1.3, 0.2],\n",
    "       [4.9, 3.6, 1.4, 0.1],\n",
    "       [4.4, 3. , 1.3, 0.2],\n",
    "       [5.1, 3.4, 1.5, 0.2],\n",
    "       [5. , 3.5, 1.3, 0.3],\n",
    "       [4.5, 2.3, 1.3, 0.3],\n",
    "       [4.4, 3.2, 1.3, 0.2],\n",
    "       [5. , 3.5, 1.6, 0.6],\n",
    "       [5.1, 3.8, 1.9, 0.4],\n",
    "       [4.8, 3. , 1.4, 0.3],\n",
    "       [5.1, 3.8, 1.6, 0.2],\n",
    "       [4.6, 3.2, 1.4, 0.2],\n",
    "       [5.3, 3.7, 1.5, 0.2],\n",
    "       [5. , 3.3, 1.4, 0.2]]), 'versicolor': np.array([[7. , 3.2, 4.7, 1.4],\n",
    "       [6.4, 3.2, 4.5, 1.5],\n",
    "       [6.9, 3.1, 4.9, 1.5],\n",
    "       [5.5, 2.3, 4. , 1.3],\n",
    "       [6.5, 2.8, 4.6, 1.5],\n",
    "       [5.7, 2.8, 4.5, 1.3],\n",
    "       [6.3, 3.3, 4.7, 1.6],\n",
    "       [4.9, 2.4, 3.3, 1. ],\n",
    "       [6.6, 2.9, 4.6, 1.3],\n",
    "       [5.2, 2.7, 3.9, 1.4],\n",
    "       [5. , 2. , 3.5, 1. ],\n",
    "       [5.9, 3. , 4.2, 1.5],\n",
    "       [6. , 2.2, 4. , 1. ],\n",
    "       [6.1, 2.9, 4.7, 1.4],\n",
    "       [5.6, 2.9, 3.6, 1.3],\n",
    "       [6.7, 3.1, 4.4, 1.4],\n",
    "       [5.6, 3. , 4.5, 1.5],\n",
    "       [5.8, 2.7, 4.1, 1. ],\n",
    "       [6.2, 2.2, 4.5, 1.5],\n",
    "       [5.6, 2.5, 3.9, 1.1],\n",
    "       [5.9, 3.2, 4.8, 1.8],\n",
    "       [6.1, 2.8, 4. , 1.3],\n",
    "       [6.3, 2.5, 4.9, 1.5],\n",
    "       [6.1, 2.8, 4.7, 1.2],\n",
    "       [6.4, 2.9, 4.3, 1.3],\n",
    "       [6.6, 3. , 4.4, 1.4],\n",
    "       [6.8, 2.8, 4.8, 1.4],\n",
    "       [6.7, 3. , 5. , 1.7],\n",
    "       [6. , 2.9, 4.5, 1.5],\n",
    "       [5.7, 2.6, 3.5, 1. ],\n",
    "       [5.5, 2.4, 3.8, 1.1],\n",
    "       [5.5, 2.4, 3.7, 1. ],\n",
    "       [5.8, 2.7, 3.9, 1.2],\n",
    "       [6. , 2.7, 5.1, 1.6],\n",
    "       [5.4, 3. , 4.5, 1.5],\n",
    "       [6. , 3.4, 4.5, 1.6],\n",
    "       [6.7, 3.1, 4.7, 1.5],\n",
    "       [6.3, 2.3, 4.4, 1.3],\n",
    "       [5.6, 3. , 4.1, 1.3],\n",
    "       [5.5, 2.5, 4. , 1.3],\n",
    "       [5.5, 2.6, 4.4, 1.2],\n",
    "       [6.1, 3. , 4.6, 1.4],\n",
    "       [5.8, 2.6, 4. , 1.2],\n",
    "       [5. , 2.3, 3.3, 1. ],\n",
    "       [5.6, 2.7, 4.2, 1.3],\n",
    "       [5.7, 3. , 4.2, 1.2],\n",
    "       [5.7, 2.9, 4.2, 1.3],\n",
    "       [6.2, 2.9, 4.3, 1.3],\n",
    "       [5.1, 2.5, 3. , 1.1],\n",
    "       [5.7, 2.8, 4.1, 1.3]])}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 14.1 Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([ 1, -1,  1])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf2pm1 = lambda b:2*b-1 #convert T/F to 1 -1\n",
    "b=True\n",
    "tf2pm1(b)\n",
    "b = np.array([True,False,True])\n",
    "tf2pm1(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "numTP = lambda y,yhat: sum([1 for i in range(len(y)) if y[i] == True and yhat[i] == True])\n",
    "numFN = lambda y,yhat: sum([1 for i in range(len(y)) if y[i] == True and yhat[i] == False])\n",
    "numFP = lambda y,yhat: sum([1 for i in range(len(y)) if y[i] == False and yhat[i] == True])\n",
    "numTN = lambda y,yhat: sum([1 for i in range(len(y)) if y[i] == False and yhat[i] == False]) \n",
    "confusion_matrix = lambda y,yhat: np.vstack([[numTP(y,yhat),numFN(y,yhat)],[numFP(y,yhat),numTN(y,yhat)]])\n",
    "error_rate = lambda y,yhat: (numFN(y,yhat) + numFP(y,yhat)) / len(y)\n",
    "error_rate2 = lambda y,yhat: np.average(y != yhat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "# y,yhat = np.random.choice(a = [True, False], size = (100)),np.random.choice(a = [True, False], size = (100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pseudo-seed to test between skl and Julia\n",
    "y = np.array([False,  True, False,  True, False,  True, False,  True,  True,\n",
    "        True,  True,  True, False,  True,  True,  True,  True, False,\n",
    "       False,  True,  True,  True,  True,  True,  True, False, False,\n",
    "       False,  True,  True,  True,  True,  True, False,  True,  True,\n",
    "        True,  True,  True,  True,  True,  True,  True,  True, False,\n",
    "        True,  True, False, False,  True,  True,  True,  True, False,\n",
    "        True, False,  True, False,  True, False,  True,  True,  True,\n",
    "        True, False,  True,  True, False, False, False, False,  True,\n",
    "        True,  True, False, False,  True, False,  True,  True, False,\n",
    "        True, False,  True, False, False, False,  True, False, False,\n",
    "        True,  True, False, False,  True,  True, False, False,  True,\n",
    "       False])\n",
    "yhat = np.array([False,  True,  True,  True,  True,  True, False, False, False,\n",
    "       False, False, False, False,  True,  True,  True, False, False,\n",
    "        True,  True, False,  True,  True, False,  True, False,  True,\n",
    "       False,  True,  True, False,  True, False, False,  True, False,\n",
    "       False, False, False,  True,  True, False, False,  True, False,\n",
    "        True,  True, False, False,  True, False, False,  True,  True,\n",
    "        True,  True,  True, False, False, False,  True, False,  True,\n",
    "       False, False, False, False, False, False,  True, False,  True,\n",
    "       False,  True, False,  True, False, False, False, False, False,\n",
    "       False, False,  True,  True, False,  True, False,  True,  True,\n",
    "       False,  True,  True, False,  True,  True, False,  True,  True,\n",
    "        True])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[32, 30],\n",
       "       [15, 23]])"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[32, 30],\n",
       "       [15, 23]])"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.45"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.45"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "confusion_matrix(y,yhat)\n",
    "np.flip(sklm.confusion_matrix(y,yhat)) #sklm implementation is flipped relative to VMLS\n",
    "error_rate(y,yhat)\n",
    "error_rate2(y,yhat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x129fa3cc0>"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFyCAYAAAA+gYtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VOXd//HPhCRmYyAUCLJHdloNYAxKZa+4VMVQAQEF6wJUJ6AiCAFcIKA/hYYnxGJEFMSILBYMbcUIuDxtCYsasWxBQMAIAymBIRAIYeb3B4+p00ASnDCTc8775TV/cJ9z7ty5Lq6PX773OWdsHo/HIwCAIQQFegEAgKojtAHAQAhtADAQQhsADITQBgADIbQBwEAIbQAwEEIbAAyE0AYAAyG0AcBACG0AMBBCGwAMJDiQP/xcwd5A/njUQKXfrA/0ElBDhfd+xKfrqyNvQupf4/McvqLSBgADCWilDQB+4z4f6BVUC0IbgDV43IFeQbUgtAFYg9scoU1PGwAMhEobgCV4aI8AgIGYpD1CaAOwBpNU2vS0AcBAqLQBWAP3aQOAgZikPUJoA7AGNiIBwDjMcssfG5EAYCBU2gCsgfYIABiISdojhDYAa+CWPwAwEJNU2mxEAoCBUGkDsAY2IgHAQEzSHiG0AViDSSptetoAYCBU2gAswePhlj8AMA562gBgICbpaRPaAKzBJJU2G5EAYCBU2gCsgXePAICBmKQ9QmgDsAY2IgHAQExSabMRCQBXgNPp1FNPPaWuXbuqc+fOGjlypHbv3l12fMeOHXrggQfUqVMn9erVSwsWLKjSvIQ2AGtwu33/VJHH49Gjjz6qw4cPa8GCBVqxYoXCwsL04IMP6tSpUzp27JgefPBBtWjRQu+//77Gjh2rtLQ0LVu2rNK5aY8AsAY/9rQLCgrUqlUrjRkzRrGxsZKkxx57TP3791deXp42btyokJAQPf/88woODlarVq20f/9+vf766xo0aFCFc1NpA7AEj+e8z5+qatCggVJTU8sCu6CgQAsWLFDDhg3Vtm1bbdmyRfHx8QoO/k/d3LVrVx08eFBOp7PCuam0AeAKmjhxolauXKnQ0FDNmzdPkZGRcjqdat26tdd5DRs2lCQdOnRIMTExl5yP0AZgDdXQHnG5XHK5XOXG7Xa77Hb7Ra95+OGHNWzYML377rt6/PHHlZmZqTNnzig0NNTrvB//fPbs2QrXQGgDsIZquOVv0aJFSk9PLzfucDiUlJR00WvatGkjSZoxY4a+/vprLV68WGFhYSopKfE678c/R0REVLgGQhuANVRDpT1ixAglJiaWG//vKvvIkSPauHGj7rzzTtlsNklSUFCQWrduLafTqUaNGunIkSPlrpGkRo0aVbgGQhuANVRDpV1RG+SnDh06pKefflpXX3214uPjJUnnzp3T9u3b1bNnT8XExCgzM1OlpaVlm5E5OTlq2bKlGjRoUOHc3D0CANXs2muvVdeuXfXss89qy5YtysvL0zPPPKPjx4/rwQcf1O9+9zsVFxcrOTlZ3377rVatWqWFCxdq1KhRlc5NpQ3AGvx4n3ZQUJDmzp2rWbNm6YknntDJkycVHx+vzMxMNWvWTJK0YMECzZgxQ4mJiWrQoIHGjRunAQMGVDq3zePxeK70L3Ap5wr2BupHo4Yq/WZ9oJeAGiq89yM+XV/8UfkNxMtew60On+fwFZU2AGswyVv+6GkDgIFQaQOwBpNU2oQ2AGswyfu0CW0A1kClDQAGYpJKm41IADAQKm0A1kB7BAAMxCTtEUIbgDVQaQOAgZgktNmIBAADodIGYA2BezdetSK0AViDSdojhDYAazBJaNPTBgADodIGYA3cpw0ABmKS9gihDcAauHsEAAzEJJU2G5EAYCBU2gCswSSVNqENwBq4ewQAjMPjZiMSAIzDJO0RNiIBwECotAFYAz1tADAQetoAYCAm6WkT2n60bNXf9M7yVcr/wamrYxpo8IA7df/A/rLZbJKkrdt2au78t7Ujb48kqWO71nryD79Xh7atA7lsXGErPs/Vu+u/VP6/T6hRvdoa1KOzhvbpUvb34pjrlP74/qf6x7Z9OnuuVDe0a66nB/ZWswbRAV45AoGNSD9ZvHSlps9KV++bb1TaS8/qjn699crc+cpYuESStHP3Xj3omCC3263nnxmj5yaM0eniMxo26intyPs2wKvHlfLOui2aseRj9YxrpdQ/3KPbb+ig2Ss+0fy/bZAknXe79djcFdqcd1Dj7u2tqfffqv3OQj36x6UqKj4b4NUbjNvt+6cGoNL2A7fbrTfeWa7f9uutJ//wkCTp112v18H8Q8pc8YFG/36o3np3hX4RHa15s6YpNDT0wjkJXdTv3gf19tJVenHq04H8FXAFuN0evfXRRt2R0FFjE3tKkrp1jNX3R49rySdfauRvu+njL3Zp58EjWpI8XB2ax0iSurRuqruffUPLP8/V72/tGshfwVh4YRSqymazaf6cGYqKjPAavyo0RCUl5yRJba5pqV91aFsW2JIUERGuq2Ma6GjBv/26XviHzSa9NnaQIsNCvcZDQ4JVUnpekvSPbfvUpH6dssCWpJjo2urUqon+95u9hPblqCGVsq8qDe2SkhKtWbNGW7Zs0aFDh3T27FlFRESoUaNGSkhIUL9+/RQcTPZXxGazqW2rWEmSx+PRCddJrf3sn8pas04PDEqUJD3ywKBy1x3MP6Tde/dryIA7/bpe+IfNZlObJg0k/d/fi1NntD53t/6Ss033942XJO07fEwtY+qVu7Z5w7pa+2WeX9dreFa4e+TAgQN6+OGHVVBQoI4dO6phw4aqV6+eSkpKtHv3bn3wwQeaO3eu5s+fr6ZNm/przYa2+cutemjMRElSx3ZtNPy+xIued+bsWSWnzNJVoSEaPvji58A8tuQd1KOpSyVJHZvH6P6+10uSiorPqkn9OuXOjwy7SkVnSvy6RtQMFYb2Cy+8oNjYWK1cuVJRUVHljhcVFenJJ5/U9OnTlZGRccUWaSbNmzXWW+n/T4ePFOhPC97R4IfH6L03/kf16/3nTgDXySKNmfiCvtmepzkzp+jqRg0DuGL4Q/OG0XrjqfvkLHTptb/8U0NfXKzMSQ/I7fHIdolrLjWOS7DCwzVffPGFli1bdtHAlqSoqCiNGzdOQ4cOvSKLM6NGDRuoUcML/yS+rmM7/fa+R/T+6jUaNWKIJOnA9z/o8QnP69DhI0pNmaxev6ZnaQUx0bUVE11bkvSrlo3V/7k39Oe/b1Xt8ItX1KfOnFVU+FX+XqaxmaQ9UuEtf3a7XU6ns8IJ8vPzFRERUeE5Vuc6WaTVH63X9z8c9hpv0ayJoiIjdNhZIEn6ZvsuDRv1pI6fOKE30l5U7+43BmK58BPX6TP668Ztyi847jXeIiZaUWGhchaeVMtG9XTgSGG5aw8cKVRso1/4a6mm4HG7ff7UBBWG9r333quJEyfqvffe0969e3X69GmVlpbq9OnT+u6777Rs2TJNnjxZAwYM8Nd6DWvqjD/q7ff+7DX21TfbVXTqtNq3uUb79n+vkU9OVmREhDIzUtXpVx0CtFL403OL1mjx2i1eY7l78lV0pkTtmjVUt46xOnCkUDsP/qd4chaeVO6eH/TrX8b6e7nG5vb4/qkBKmyPJCUlyWaz6eWXX1ZxcXG545GRkRo2bJjGjh17xRZoBvbaURox5Hd6690VioiI0A2dr9W+A9/r9YVL1KFtK91zxy0a9dRknTpdrGfGjFLBv4+p4N/Hyq6Pioosu/sE5mGPCNPwW27QwuxNiggLVXzbZvru8DHN/zBH7Zs11N03/UpBQTa99dFGJaW/r8f7d1dYSLBe+8s/Va92hAb27BToXwEBYPN4Kr/jvKSkRDt37pTT6VRxcbHCwsLUqFEjtW/f3uu+4st1rmDvz77WaNxut5au/KuWrvyrDuT/oLp2u27pfbOSHn1AZ0vOqeedQy55befrOmrxvNl+XG3glH6zPtBL8Cu326Nln3+l5Z9/rYNHClUnMly/6dJWj999c1nP+sjxIr2ybL02bN8n2WyKb9NM4wb2stxj7OG9H/Hp+lMp9/u8hsgp7/g8h6+qFNpXipVCG1VjtdBG1fkc2tOG+byGyGczfZ7DVzwVA8AaashGoq94YRQAGAiVNgBrqCF3f/iK0AZgDVZ4IhIATMPPlXZRUZHS0tK0du1aFRYWKjY2Vo8//rj69u2riRMnauXKlRe9LikpSQ6H45LzEtoALMHfTzROmjRJu3btUkpKipo0aaIPP/xQDodDb775piZPnqxx48Z5nT937lytXbtWAwcOrHBeQhsAqtnRo0eVnZ2tjIwMdevWTZI0evRobdiwQStWrNBNN92k2rVrl52/ZcsWLV++XK+99ppiYmIuNa0kQhuAVfixPRIeHq758+erS5cuXuM2m00nTpzwXpbbrenTp+uWW25Rz549K52b0AZgDdUQ2i6XSy6Xq9y43W6X3W4v+3NUVJR69OjhdU5ubq5ycnI0ZcoUr/E1a9YoLy9Pc+bMqdIaCG0A1lANd48sWrRI6enp5cYdDoeSkpIued2ePXvkcDgUFxenwYMHl5vz1ltvVWxs1d4vRGgDsIZqqLRHjBihxMTy3yT10yr7v23evFkOh0ONGzdWRkaGQkJCyo4dOHBAubm5euKJJ6q8BkIbAKrov9sglcnKylJycrISEhKUlpZW7gtl1q5dqwYNGighIaHKcxLaACzB4+f7tFevXq0JEyborrvu0syZM70q7B9t2bJFCQkJqlWrVpXnJbQBWIMfQ/vw4cOaOnWqunbtqvHjx+v48f98O1FISIjq1q0rSdqxY0e5HndlCG0A1uDHh2uys7NVXFysnJwcde/e3etYly5dtGTJEkkX7ueuU6fOZc1NaANANRs+fLiGDx9e6Xn/+te/LntuQhuANfCWPwAwEEIbAIwjgN+sWK0IbQDWYJJKm68bAwADodIGYA0mqbQJbQCW4O8nIq8UQhuANRDaAGAg5vheXzYiAcBIqLQBWAI9bQAwEkIbAAyEnjYAwN+otAFYAj1tADASk7RHCG0AlkClDQBGYpJKm41IADAQKm0AluAxSaVNaAOwBkIbAIyDShsAjMQkoc1GJAAYCJU2AEugPQIABkJoA4CBmCW06WkDgIFQaQOwBo8t0CuoFoQ2AEswS3uE0AZgCR43lTYAGIZZKm02IgHAQKi0AViCh41IADAOs7RHCG0AlmCWjUh62gBgIFTaACzBY47v9SW0AViDWdojhDYASyC0AcBAzNIeYSMSAAyEShuAJdAeAQAD4YlIADAQszwRSU8bgCW4PTafP5ejqKhIM2fOVJ8+fdS5c2cNGDBA69atu+i58+bNU7t27VRaWlrpvIQ2AFwBkyZN0qeffqqUlBStWrVK/fr1k8Ph0IYNG7zO27p1q9LT06s8L6ENwBI8HpvPn6o6evSosrOzlZycrG7duqlFixYaPXq0EhIStGLFirLzTp8+rfHjxys+Pr7KcxPaACzB47b5/Kmq8PBwzZ8/v1wY22w2nThxouzPM2bMUNu2bdW/f/8qz01oA7AEj8f3T1VFRUWpR48eioqKKhvLzc1VTk6OevXqJUn6+OOP9dlnn2natGmX9Xtw9wgAVJHL5ZLL5So3brfbZbfbL3ndnj175HA4FBcXp8GDB8vpdGrq1Kl6+eWXFR0dfVlrILQBWEJ1PFyzaNGii24aOhwOJSUlXfSazZs3y+FwqHHjxsrIyFBwcLAmTpyo22+/XT169LjsNdg8nsA9kX+uYG+gfjRqqNJv1gd6Caihwns/4tP1/7rmTp/X0Dz33cuqtLOyspScnKyEhASlpaUpKipK+fn56tOnj8LCwhQUdKFDXVpaqpKSEkVERGjUqFEaPXr0JddApQ3AEqrjicjK2iA/tXr1ak2YMEF33XWXZs6cqZCQEElSTEyMsrOzvc7Nzs7WrFmz9P7776tevXoVzktoA7AEf/YUDh8+rKlTp6pr164aP368jh8/XnYsJCRELVq08Dr/F7/4hSSpefPmCg6uOJYJbQCoZtnZ2SouLlZOTo66d+/udaxLly5asmTJz56bnjZqFHrauBRfe9q5Le72eQ2d9mf5PIevqLQBWAJv+QMAAzHLN9cQ2gAs4XLf0ldTBTS0wxt3r/wkWMra6G6BXgJqqF5O33raZkGlDcAS6GkDgIHQHgEAAzHJPiSvZgUAI6HSBmAJtEcAwEDYiAQAA3EHegHVhNAGYAkemaPSZiMSAAyEShuAJbhNcs8foQ3AEtwmaY8Q2gAswSw9bUIbgCWY5e4RNiIBwECotAFYAu0RADAQs7RHCG0AlmCW0KanDQAGQqUNwBLoaQOAgbjNkdmENgBr4IlIADAQk7x6hI1IADASKm0AlmCWW/4IbQCW4LbR0wYAwzBLT5vQBmAJZmmPsBEJAAZCpQ3AEni4BgAMhIdrAMBAzLIRSU8bAAyEShuAJdDTBgADMcstf4Q2AEswS0+b0AZgCWZpj7ARCQAGQqUNwBLoaQOAgRDaAGAgHpP0tAltAJZglkqbjUgAMBBCG4AluKvhczmKioo0c+ZM9enTR507d9aAAQO0bt26cuedPXtW/fv31/Lly6s0L6ENwBI81fC5HJMmTdKnn36qlJQUrVq1Sv369ZPD4dCGDRvKzjl58qQee+wx7dy5s8rzEtoALMFt8/1TVUePHlV2draSk5PVrVs3tWjRQqNHj1ZCQoJWrFghSVq/fr3uvvtuFRYWXtbvQWgDQDULDw/X/PnzFR8f7zVus9l04sQJSdInn3yiYcOG6b333rusubl7BIAlVMfdIy6XSy6Xq9y43W6X3W4v+3NUVJR69OjhdU5ubq5ycnI0ZcoUSdL06dN/1hoIbQCWUB2hvWjRIqWnp5cbdzgcSkpKuuR1e/bskcPhUFxcnAYPHuzTGghtAJZQHW/5GzFihBITE8uN/7TK/m+bN2+Ww+FQ48aNlZGRoZCQEJ/WQGgDsITqeMvff7dBKpOVlaXk5GQlJCQoLS1NUVFRPq+BjUgAuAJWr16tCRMm6Pbbb1dGRka1BLZEpQ3AIvz5GPvhw4c1depUde3aVePHj9fx48fLjoWEhKhu3bo/e25CG4Al+POba7Kzs1VcXKycnBx1797d61iXLl20ZMmSnz03oQ3AEtx+jO3hw4dr+PDhVT5/165dVT6XnjYAGAiVNgBLMMurWQltAJbAt7EDgIFQaQOAgVTHwzU1ARuRAGAgVNoALMGft/xdSYQ2AEswR2QT2gAsgo1IADAQs7RH2IgEAAOh0gZgCeaoswltABZBTxsADISeNgDA76i0AViCOepsQhuARdDTBgAD8Zik1ia0AViCWSptNiIDoHHjRjrq3Ka+fby/8PN/5qSotCS/3GfC+McDtFJccUFBavLI7brhs9nqvm+xum6cq1YvjFCtyLCyU+rc1FGds6br5j1v66atr6t1yu9VKyo8gItGIFFp+1nTpo31t79mKjq6brljneJ+qTVr1uvFl9K8xr/b/72/lgc/i510n5qNvksH/5Sl4xu2K6JNE7UcN1D2+Lb66s4pqt2pleKWTtGxT3K1fWSqrmoUrdjkoYps30xf3zst0Ms3FLPc8kdo+4nNZtP999+rl1+aqqCgi/8D59prO+jFl9L0j39u9vPqEAhB4aFqNvoufZ/xF+17cYkkqfDTr3Wu4IQ6vvaE6v76l2o8op/OHj6mbQ/Plqf0vCTJ4/ao/ZzHFNGuqU7v4n/oVWWOyCa0/ea66zpq3qsvad68RVr/yd+1Omux1/FWrVrKbq+trVu3B2iF8LfgOlE6/N4nOpK1wWv81K6DkqSrYqK15/m3FVwnoiywJcl99pwkKeiqUP8t1gTMUmnT0/aTAwfy1a7DzRr/zDSdPl1c7nhc3C8lSYmJd2jvt5tUfOo7bd70kW67tbe/lwo/KTl8THnjX1fR1r1e4/VvvUGSVLTjgM7mF+jU9gOSpFoRYYruca2umTREJzbtLHcdKuauhk9NQGj7SWHhceXnH7rk8U7/F9r169fToyPHaeCgR3Xs34X6YNUi3dqvl59WiUCzX99GzcckqmDNZp3avr9s3FYrSDfvXqi45c+qVmSYvp26MHCLREDRHqkh3lr4nnJyvtCHa9bL47nwz7g1H32ir75cq+efH6+Psj8N7AJxxdXt1lG/WjhBZw44tXPsn7wP1grS1iEzZQsNVrNRv1XnD6Zp69CZOv6PbYFZrAFZ5j7toUOHymar2tcYZ2Zm+rwgq9q374D27TvgNVZaWqrs7M/02B9GBGhV8JdGg3up7SsjdWrnQW0dOkOlx4u8jntKSlX4+VZJUuHnW5Xw9zlqPnYAoX0Zakp7w1eVhnbPnj01Z84cXXPNNbruuuv8sSZLSky8Q+fPn1dW1kde4+HhYSooOBagVcEfWo4fpJZPD9S/132l7Y/8UedPnyk7Vv+OBJUcPSHX5l1lY56SUp3acUARrRoHYrmGZZlKe9SoUYqKitLs2bOVkZGhpk2b+mNdljPkvnvU/eYb9ckn/9DJkxeqrMjICP32jr769LN/Bnh1uFKaOfqr5dMDdShznXY9/brk9q4HW4wdIFtwLW255ZmyY8H2CNmvb0OVbVFV2ogcNmyYEhISNGfOnCu9Hst66aW5ql07Un/7S6b6979NAwferfXr3ldkZISee/6VQC8PV0B4bCPFTrxPp/K+1+Gln6pOfFvVSWhf9glpUEf7Xl6qyA7N9cs3nlK9vp3VcMDN6rTyBdlCgrXv5WWB/hUMxSx3j1R5I3LatGnato3/s18pX371jfr+ZqCmvTBBb7w+W7Vq1dLn/5uj3z80Vnv37q98AhhO/TsSFBQSrMi2TdU5a3q547uemqdDmeu19b4Zavn0QHXMeEKe824d//u/tG1kqor3/BCAVRuX22OO9ojN4wncbxIc2iRQPxo11NroboFeAmqoXs7lPl1/f4sBPq/hnf1/9nkOX3HLHwBL4IlIAIDfUWkDsATL3PIHAGZQU+7+8BWhDcASzNLTJrQBWIJZ2iNsRAKAgVBpA7AEetoAYCABfI6wWhHaACzBLBuR9LQBwECotAFYAj1tADAQbvkDAANxy+Pz5+fKyMjQkCFDvMb27t2rRx99VPHx8br55ps1ZcoUuVyuSucitAFYgsfj8fnzc2RmZio1NdVrrKSkRCNHjlRwcLCWLl2qOXPmaPPmzUpOTq50PtojAHAFOJ1OPffcc9q4caNiY2O9ju3atUsHDx7Uq6++qlatWkm68A1hs2fPrnReKm0AluDvrxvbtm2bIiMjlZWVpbi4OK9j0dHRstlsWr58uUpKSnTs2DF9/PHH6tSpU6XzUmkDsAR/b0T26dNHffr0ueixpk2basqUKZo9e7YyMzPldrvVpk0bvf3225XOS6UNwBKqYyPS5XLp+++/L/epygbiT5WUlGjHjh36zW9+o6VLl2r+/PnyeDwaO3asSktLK7yWShuAJVTHY+yLFi1Senp6uXGHw6GkpKQqz7Nw4UJt2rRJH374oYKDL8RwixYt1K9fP61du1a33XbbJa8ltAGgikaMGKHExMRy43a7/bLm+eKLL9ShQ4eywJYuhHZ0dLT2799f4bWENgBLqI53j9jt9ssO6IuJiYnRpk2b5Ha7FRR0oUvtdDp1/PhxtWzZssJr6WkDsARPNfxXXYYNG6b8/HxNnTpVe/bsUW5ursaMGaM2bdqod+/eFV5LaAOwBLfH4/OnurRr106LFy/WwYMHdd9998nhcCg2NlYLFy5UaGhohdfaPAF8yWxwaJNA/WjUUGujuwV6CaihejmX+3R9jyZ9fV7D5/nrfJ7DV/S0AViCOV4XRWgDsAizfAkCoQ3AEghtADAQs3xHJHePAICBUGkDsATaIwBgIGb5ujFCG4AlmKWnTWgDsASztEfYiAQAA6HSBmAJtEcAwEDM0h4htAFYglnuHqGnDQAGQqUNwBKq833YgURoA7AEs7RHCG0AlkClDQAGYpZKm41IADAQKm0AlkB7BAAMxCztEUIbgCVQaQOAgZil0mYjEgAMhEobgCV4PO5AL6FaENoALIG3/AGAgZjlfdr0tAHAQKi0AVgC7REAMBCztEcIbQCWwMM1AGAgPFwDAPA7Km0AlkBPGwAMhLtHAMBAzFJp09MGAAOh0gZgCdzyBwAGYpb2CKENwBLYiAQAAzFLpc1GJAAYCJU2AEtgIxIADMQs7x4htAFYApU2ABgIG5EAgCrJyMjQkCFDvMZmz56tdu3alfuUlpZWOBeVNgBLCFRPOzMzU6mpqercubPX+K5duzRo0CCNGTPGazw4uOJYJrQBWIK/2yNOp1PPPfecNm7cqNjY2HLH8/Ly1Lt3bzVo0OCy5qU9AsASPB6Pz5/LsW3bNkVGRiorK0txcXFex1wulw4dOqTWrVtf9u9BpQ0AVeRyueRyucqN2+122e12r7E+ffqoT58+F50nLy9PkrR69WpNnjxZ586dU0JCgsaNG6eGDRtWuIaAhnZpSX4gfzwACzlXDXkzd+5cpaenlxt3OBxKSkqq8jw/hnZUVJTS0tJ09OhRpaam6oEHHtCqVasUHh5+yWttHrPcBwMAV9jlVNo/NXHiRO3fv19LliyRdKFV43K5VKdOnbJznE6nevbsqVdeeUV33XXXJeeiPQIAVVRZOFeVzWbzCmxJiomJUd26dXXo0KEKr2UjEgD8LCUlRffcc4/X2MGDB1VYWFjp5iShDQB+dtttt2n37t1KSUnRd999p02bNsnhcOi6665Tr169KryW9ggA+Fl8fLxee+01paenKzExUaGhoerbt6/Gjx+voKCKa2k2IgHAQGiPAICBENoAYCCENgAYCKEdQG63W2lpaerevbvi4uL00EMPaf/+/YFeFmqQi73SE9ZGaAfQq6++qiVLliglJUVLly5VrVq19PDDD+vs2bOBXhpqgB9f6Qn8FKEdICUlJXrzzTflcDjUs2dPtW/fXqmpqSooKNCHH34Y6OUhgJxOp0aPHq1Zs2Zd9JWesDZCO0B27Nih06dP68Ybbywbi4qKUseOHbVly5YArgyBVtErPQEergkQp9Mp6cL7Bn5ulWjdAAAA8UlEQVSqYcOGlb57AOZW0Ss9ASrtACkuLpYkhYaGeo2HhoaqpKQkEEsCYACEdoCEhYVJUrmALikpUURERCCWBMAACO0AufrqqyVJR44c8Ro/cuRIuZYJAPyI0A6Q9u3bKyoqSps2bSobKyoq0vbt25WQkBDAlQGoydiIDJDQ0FDdf//9Sk1NVf369dW0aVPNnj1bMTEx6tevX6CXB6CGIrQDaMyYMTp//ryeffZZFRcX6/rrr9cbb7xRbnMSAH7Eq1kBwEDoaQOAgRDaAGAghDYAGAihDQAGQmgDgIEQ2gBgIIQ2ABgIoQ0ABkJoA4CB/H/2PKh63hyyPQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#optional confusion matrix vis template\n",
    "import seaborn as sn\n",
    "import pandas as pd\n",
    "\n",
    "df_cm = pd.DataFrame(confusion_matrix(y,yhat) )\n",
    "sn.heatmap(df_cm, annot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 14.2 Least Squares Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[46,  4],\n",
       "       [ 7, 93]])"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.07333333333333333"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.07333333333333333"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#fHat can be made via a regression model, and the > comparator\n",
    "fTilde = lambda x: x@beta + v\n",
    "fHat = lambda x: fTilde(x) > 0\n",
    "\n",
    "\n",
    "iris = np.vstack([irisD[\"setosa\"],irisD[\"versicolor\"],irisD[\"virginica\"]])\n",
    "#if k == virginica: y[k] = True (1) ; else False (0)\n",
    "y = np.hstack([np.full(50, False),np.full(50, False),np.full(50, True)])\n",
    "A = np.hstack([np.ones((150,1)), iris])\n",
    "theta = npl.lstsq(A,2*y-1)[0]\n",
    "yhat = np.matmul(A,theta) > 0 #regression classifier \n",
    "\n",
    "confusion_matrix(y,yhat)\n",
    "(C[0,1] + C[1,0]) / len(y)\n",
    "np.average(y != yhat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 14.3 Multi-class Classifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 306,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 4\n",
    "y = np.random.randint(K,size=100)\n",
    "yhat = np.random.randint(K,size=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 307,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6., 11.,  8.,  9.],\n",
       "       [ 8.,  7.,  5.,  6.],\n",
       "       [ 4.,  3.,  6.,  3.],\n",
       "       [ 7.,  7.,  3.,  7.]])"
      ]
     },
     "execution_count": 307,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.74"
      ]
     },
     "execution_count": 307,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.74"
      ]
     },
     "execution_count": 307,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "error_rate = lambda y,yhat: np.average(y != yhat)\n",
    "def confusion_matrix(y, yhat, K):\n",
    "    C = np.zeros((K,K))\n",
    "    for i in range(K):\n",
    "        for j in range(K):\n",
    "            C[i,j] = sum(np.logical_and(y == i, yhat == j))\n",
    "    return C\n",
    "C = confusion_matrix(y,yhat,K)\n",
    "C\n",
    "error_rate(y,yhat)\n",
    "1-np.sum(np.diag(C))/np.sum(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 308,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x12b3310b8>"
      ]
     },
     "execution_count": 308,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAFyCAYAAAD78xH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8jOf6P/DPZJmsJrYs1SDEVmoJkqAcSoO26ElPW0sRsRRtQlpqrdiporRS5SgaO9VaWorztbRah4i0tkRCgwjZkGayzmSZ3x9+0jOdiGVi7uS5P+++5vVq7nlm5jPpuFy9nmVUBoPBACIiUiwr0QGIiOjpYqEnIlI4FnoiIoVjoSciUjgWeiIihWOhJyJSOBZ6IiKFY6EnIlI4FnoiIoVjoSciUjgWeiIihWOhJyJSOBuRL76q7mCRL19pjPhtjugIlYY2OFh0hErDoXND0REqFcfJ6816fOHtRLMz2Naumv9N2NETESmc0I6eiMhiSopFJxCGhZ6I5GAoEZ1AGBZ6IpJDibyFnjN6IiKFY0dPRFIwcHRDRKRwEo9uWOiJSA4Sd/Sc0RMRKRw7eiKSA4+jJyJSOIlHNyz0RCQH7owlIlI2mQ+v5M5YIiKFY6EnIjmUlJh/e0KrV6/GwIEDjdaSk5MxevRotG3bFp06dcLixYtRVFRU7vP0798fTZs2Nbr9/XnLwtENEclB0Ohm8+bNWLZsGXx8fErX9Ho9RowYgQYNGmDbtm24ceMGpk2bBhsbG7z//vtlPo/BYEBCQgIWLlyILl26lK7b2to+NAMLPRHJwcKHV6alpWHmzJk4deoUGjRoYHTfwYMHcfPmTezYsQMuLi5o0qQJJk6ciAULFmDs2LGwt7c3eb7k5GTk5eWhdevWcHV1fawsHN0QkRwMJebfHsPFixfh5OSEvXv3onXr1kb3RUdH47nnnoOLi0vpmr+/P/Ly8nDx4sUyny8+Ph62traoX7/+Y791dvRERE9B9+7d0b179zLvS0tLg4eHh9Gam5sbACA1NbXMx8THx8PFxQVTp07FqVOn4OzsjF69emHs2LFQq9XlZmGhJyI5VMBx9FqtFlqt1mRdo9FAo9E88vMUFBTAycnJaO1+sdbpdGU+5vLly8jLy4Ofnx9GjRqF2NhYLFq0CMnJyVi8eHG5r8dCT0RyqICdsZGRkYiIiDBZDwkJQWho6CM/j729PfR6vdHa/Z8dHR3LfMySJUtQUFAAZ2dnAECTJk1gY2ODCRMmYOLEiXB3d3/g67HQE5EcKqCjDwoKQmBgoMn643TzAODh4YG4uDijtfT09NL7ymJjY1Na5O9r2rQpACAlJYWFnoioIjzuiOZBfH198d1330Gr1ZY+36lTp+Dk5ITmzZuX+ZjAwEC0bdsWM2bMKF07d+4cbG1t4eXlVe7r8agbIpKCwVBs9q2ivPTSS3B3d8f777+PS5cu4ciRI1i6dCmCg4NLZ/W5ubnIyMgofczLL7+M7du345tvvsGNGzewb98+fPLJJxgyZAiqV69e7uuxoyciOVSia93Y2dnhq6++wpw5c/DWW29Bo9Ggf//+eO+990q3WbduHSIiIhAfHw8AGDVqFOzs7LB27VrMmTMHrq6uCAoKwpgxYx76eiqDwWB4au/mIVbVHSzqpSuVEb/NER2h0tAGB4uOUGk4dG4oOkKl4jh5vVmPL4jZa3YG+7b9zH4OEdjRE5EcKlFHb2mc0RMRKRw7eiKSA79KkIhI4SQe3bDQE5Ec+FWCREQKx45ePm4+3vCf0h9ubRqiMFeHGz+dw8n5W5F/2/SCRUqVlnEbgUPGYuncaejo61PmNuOnzoWzkyPmfzTBwunEsOvZB/avvQFrNw8UZ6RD9+MeFPzwHSDuKGRBVLDx6QYbn+5QubjCkH0HRTFHUBRzWHQwegJSHnVTu6UX+u2YjpLCIvxnzApELdqBuv9oid7rPhAdzWJS0jIwKmwatNk5Zd5fXFyMBZ+uxOGfT1g4mTh2vfrCOfRDFF04i+wFM6A/8RMcR4bA/l8P/6o2pbHtMRDqnkNRcusP6PasRFHMEdi+8BpsX+wvOtqTE/hVgqJJ2dF3mDYQdxOSsT9oCQzF9/7j6bJy8cLsIXDxckfWtTTBCZ+ekpISfH/gMJZ88RVKSsruUi9dTsTCZV8iNv4y7B5ynWslsQt4BYWx55D75TIAQOFvp2FdxxP2rwSiYOcWweksyMEZNm17oOjccegPfA0AKAFQknUbdq+PQ9HZn2C4W/Y10yu1KlyozSVdR29X3Rl1Oj6Hi5H/V1rkAeDqgWhs8h+v6CIPAAlXrmL24hXo1/slLJwxscxtpsz+BACw7avPUN2lmiXjCaVS28GQm2u0Zsj6E1bV5PkdAIBVTQ+orKxQfOV3o/WSpDiorKxg3bCloGTmqUzXurG0h3b0er0eBw4cQHR0NFJSUqDT6eDo6AgPDw/4+fmhZ8+esLGpOv9jUOu5urCytkJ+Rha6Lx8Dr17toFKpcPXgGfwSHgl9Vp7oiE/VMx5u2L99LTzcXBEVc67MbRbNmoymjRqUeZ+SFezdCaeQCVB3C0Bh1AnYNG0Ou+69oTtySHQ0izLkZQMAVC61jdZV1d3+//rjfV8piVduhU5KSsKIESNw+/ZtNG/eHG5ubqhZsyb0ej0uX76MPXv2YMWKFVizZg08PT0tldksDrXuXRK06+KRuHHsHA6OXA6Nlzv8p7yFV70mYdc/Zyt6x5uLphpcNOV3qDIWeQDQHTsEmxatUG3CR6Vr+jOnkLvmc4GpLM+QmYbiGwmwfaEfDNo7KL4WC1V1V6h7D4OhqBAqtZ3oiE9G4tFNuYV+9uzZaNCgAXbt2mVywXsAyMnJwfvvv4+5c+di9erVTy1kRbKyvfeWb1+4jmMT1wAAbv56EXptLgJWhqJet1ZIOnpWZEQSpNpHC2DbvCVy169CUUIsbOo3hMOgYag2dS6y501TdAPwd7rdX0DdKwh2gSEAAIMuD4XHvoGqU18YCvUPeXQlxcMry3bmzBns2LGjzCIPAM7OzpgwYQIGDRr0VMI9DYW5+QCApCPG88cbP50HANR+vj4LvYRsmrWAup0/cr5YCt2Be1c5LLpwFsWpt6CZ9Qls/V9A4clfBKe0oDwt9LtWQG/nAFW1GjBkpgMlJbANGAIU5D788ZWRxB19uTtjNRoN0tLK3zl58+bNB37HYWWUdfXe0QLWauO/46xsrAEARQWFFs9E4lm53fv6tqK4C0brhRfv7cewqedl6UhCWT/nB5VbXUCXD8PtW0BxEazc60FlZYWS1Oui4z0ZQ4n5tyqq3EL/xhtvYMqUKdi2bRsSExORl5eHoqIi5OXl4dq1a9ixYwemT5+O119/3VJ5zZZ5+Ra0Senw7tfBaL1+wL0ThlKi4kXEIsGKk5MAALbPtzZat23R6t79qSkWzySSbce+sO1kfO11m/Y9YSjIQ3HSJUGp6EmVO7oJDQ2FSqXCJ598gvz8fJP7nZyc8Pbbb2P8+PFPLeDTcHL+VgR8GYpea8IQu/kIqns/A78P38TVg9HIOJsoOh4JUJx4Gbpfj8ExeAxUTk4oio+DdT0vOAwYiqI/LkN/4ifRES2qMPoQ7F4eDkPnf6I4OQE2zfxg06Ij9AcjAb1pLagSJB7dlFvoVSoVQkNDMXr0aFy6dAlpaWnIz8+Hvb09PDw80KxZs9LvN6xKEvefxo/DP0X7sED0WhMGXVYeYrccRdSiHaKjkUA5S+bC4a0hsOvZBw4DglCSkQ7d4QPI2/o1UFQkOp5FFZ87Dr2NGjbtAmDj1xuGu6nQ7V2F4rhToqM9uSo8ejEXv0qwEuBXCf6FXyX4F36VoDFzv0ow/0fzD5N1eHmc2c8hgnRnxhIRyabqnNJKRGQOzuiJiBRO4hk9Cz0RyYEdPRGRwknc0XNnLBGRwrGjJyI5cHRDRKRwEo9uWOiJSA7s6ImIFE7iQs+dsURECseOnojkINE3hP0dCz0RyUHi0Q0LPRHJQeJCzxk9EZHCsaMnIjnwOHoiIoWTeHTDQk9EcuBRN0RECidxR8+dsURECseOnojkIHFHz0JPRHLgUTdERMpmKOHOWCIiZZN4dMOdsURECseOnojkwBk9EZHCcUZPRKRwnNETEZFSsaMnIjmwoyciUjiDwfzbY8jNzcXcuXPRtWtXtGvXDu+++y6SkpIeuH1mZiYmTJgAPz8/+Pr6YsaMGcjNzTX3XQNgoSciWZSUmH97DGFhYTh48CDCw8OxY8cO1KxZEwMHDkRmZmaZ248bNw5JSUlYv349IiIicOLECYSHh1fEO2ehJyJJlBjMvz2iS5cu4eeff8acOXPQo0cPeHt7Y/bs2XB2dsaWLVtMto+JiUFUVBQWLlyIFi1awN/fH/PmzcO+fftw69Yts986Cz0RUQW7du0aAMDX17d0zdraGs2aNUNUVJTJ9tHR0ahVqxYaNWpUutauXTuoVCpER0ebnYeFnojkYCgx//aIXF1dAcCkG09OTsbdu3dNtk9PT4eHh4fRmlqtRo0aNZCamvoEb9aY0KNu3mxzQ+TLVxpn23wgOkKlsdLGVXSESsPv90LRESqVMZPNfIIKOGFKq9VCq9WarGs0Gmg0mtKfW7VqhUaNGmHWrFlYunQpateujU2bNuHSpUvw9PQ0eXx+fj7UarXJulqthk6nMzs3D68kIikYKuDwysjISERERJish4SEIDQ0tPRnW1tbREREYMqUKXjxxRdhY2ODbt264Y033sCFCxdMHm9vbw+9Xm+yrtfr4ejoaHZuFnoikkMFdPRBQUEIDAw0Wf/fbv6+Bg0aYPv27cjKyoJKpYJGo8H48ePh5eVlsq2HhwfS09ON1vR6PTIzM01GOk+CM3oiokek0Wjg6elpcvt7oc/JycHgwYNx/vx5uLi4QKPRICcnBydOnECXLl1MntfX1xcZGRlITEwsXbu/E7Z9+/Zm52ahJyI5WHBnrLOzM1QqFRYsWID4+HhcunQJY8aMQZ06ddCnTx8UFxcjIyMDBQUFAIDWrVujbdu2mDBhAs6dO4eoqCiEh4fjtddeg7u7u9lvnYWeiORgwePoAWDp0qVwdXXF4MGDERQUBE9PT6xfvx42NjZISUlB586dsX//fgCASqVCREQE6tati6CgIISGhqJTp06YNWtWhbx1zuiJSA4WvtaNm5sbPv/88zLv8/T0RHx8vNFarVq1Hri9udjRExEpHDt6IpIDv3iEiEjh+FWCREQKx46eiEjZKuLM2KqKO2OJiBSOHT0RyYGjGyIihWOhJyJSOB51Q0SkcBJ39NwZS0SkcOzoiUgKBok7ehZ6IpIDCz0RkcLxhCkiIlIqdvREJAeOboiIFI6FnohI2QwGFnoiImWTuKPnzlgiIoVjR09EcpC4o2ehJyIp8MxYIiKlY6EnIlI4eU+M5c5YIiKlY0dPRFLgjJ6ISOlY6OVj17MP7F97A9ZuHijOSIfuxz0o+OE7QKKz56p1fB5Nv5n3wPtvfboNtz7dZsFE4n1xfgMcqjmarIf5joQ2408BicRx8/GG/5T+cGvTEIW5Otz46RxOzt+K/Nta0dGejMQzeikLvV2vvnAOmYiC/buRd/IL2LRoBceRIYCdHQp2bhEdz2LyLiTiUuBUk/U6EwbCqXUj3N1zXEAqcdzqe8ChmiO2zFqH6xcSje7LzcwWlEqM2i290G/HdKScuoT/jFkBR7fq8Jv0Jqp7P4Nd/WaJjkePSc5CH/AKCmPPIffLZQCAwt9Ow7qOJ+xfCZSq0Bdn5yHndJzRmkuALzSdW+GPMYtR8MdNQcnEqNvcCwAQ9cOv0N7OEhtGsA7TBuJuQjL2By2BofheK6zLysULs4fAxcsdWdfSBCd8fJzRS0altkPJ3/4gG7L+hFW1aoISVQ4qezXqzR2FPw9HI/OHX0XHsbh6zb2QlZEpfZG3q+6MOh2fw0+Tviot8gBw9UA0rh6IFpjMTBzdyKVg7044hUyAulsACqNOwKZpc9h17w3dkUOiownlPrIf1B61kDAgXHQUIeo1b4B8bR5C10xGsw4toFKpcPZoDLbNWY8siebztZ6rCytrK+RnZKH78jHw6tUOKpUKVw+ewS/hkdBn5YmO+ERk7uilPI5ed+wQdEf/g2oTPkLN7fuhmbMEhXHnkbvmc9HRhFGpbeA+/FXc3XscumupouMIUbe5F2p5uuHy6TgsH74QOxZuRLMOLTB522yoHexEx7MYh1oaAEDXxSNRUlSMgyOX48TcLajXvTVe3TAJUKkEJ3xCJRVwq6Kk7OirfbQAts1bInf9KhQlxMKmfkM4DBqGalPnInveNKmOvLmvRt/OsHWrgdQvd4uOIsy/xy9HQW4Bki5eBQBcPh2Hm/FJmLpzHjq/8SKObDwgOKFlWNneKwu3L1zHsYlrAAA3f70IvTYXAStDUa9bKyQdPSsyIj2mhxb6QYMGQfWIf4Nv3rzZ7EBPm02zFlC380fOF0uhO7AXAFB04SyKU29BM+sT2Pq/gMKTvwhOaXk1XumI/Pgk5MddEx1FmISoOJO1y9GXkJuVU7qjVgaFufkAgKQjvxut3/jpPACg9vP1q2ShN1ThjtxcDy30Xbt2xfLly9GwYUO0atXKEpmeKis3DwBAUdwFo/XCi+cAADb1vKQr9Cq1DVy6tkHqyl2iowhTraYGbXv7I/7kRaQm3ipdV6lUsFHbIiezih47/gSyrt4b3VmrjcuDlY01AKCooNDimSoEC/2DjR49Gs7Ozli6dClWr14NT09PS+R6aoqTkwAAts+3RvH1v46Vtm1x7y+x4tQUIblEcmzeAFb2dsiOihUdRZiioiIMnjMSv+w4gshpq0vXfXr6wc7BDpf+e1FgOsvKvHwL2qR0ePfrgHNf/TWuqh/gAwBIiYoXFc0s7Ogf4u2338bx48exfPlyLFmy5GlneqqKEy9D9+sxOAaPgcrJCUXxcbCu5wWHAUNR9Mdl6E/8JDqixTk0qw8AKEi4ITiJOPnaPBxa+wN6jeqLvOw8xP56DnWf80K/0Ddw9vAZXDxe9UYV5jg5fysCvgxFrzVhiN18BNW9n4Hfh2/i6sFoZJxNfPgTVEYs9A83Z84cXLyojK4mZ8lcOLw1BHY9+8BhQBBKMtKhO3wAeVu/BoqKRMezOJva1QEARVk5gpOItXPRZvyZdhf/GPASegT1RvYdLY5sPIA9y78RHc3iEvefxo/DP0X7sED0WhMGXVYeYrccRdSiHaKj0RNQGQR+Nfqdvl1FvXSlcvW3GqIjVBorbaQ8EKxMfsX2oiNUKmNubDLr8RkB5tcb1/9Uzf/j558qIpICZ/RERAonc6GX8sxYIiKZsKMnIjkYquilGyoACz0RSUHm0Q0LPRFJwVDCjp6ISNFk7ui5M5aISOHY0RORFAzcGUtEpGyWHN2cOnUKQ4cOLfM+T09PHD582GT9+PHjGDlypMn6+vXr0alTJ7PysNATkRQsuTPWx8cHv/xifLnzhIQEvPPOOxg9enSZj4mPj0fjxo2xfv16o3UXFxez87DQExFVMLVaDVdX19KfCwsLsWDBAgQEBOCtt94q8zEJCQlo3Lix0eMqCnfGEpEUDAbzb09q48aNSElJwdSpUx+4TXx8PLy9vZ/8RcrBjp6IpCDqOPr8/HysXr0aQ4cOhbu7e5nbFBUVITExEQkJCQgMDERGRgaaNGmCsLCwCvlmPxZ6IpJCRRR6rVYLrdb0ayU1Gg00Gk2Zj9mzZw90Ot0Dd84CwPXr16HX61FQUIDw8HBYWVlhw4YNGDx4ML799ls0btzYrNws9EQkhYr45o3IyEhERESYrIeEhCA0NLTMx+zZswcBAQGoWbPmA5/X29sbp0+fhrOzM6ys7k3UFy9ejD59+mDDhg2YO3euWblZ6ImIHlFQUBACAwNN1h/Uzd+9exe///47xowZ89Dn/vtzWFlZoVGjRrh169YDHvHouDOWiKRgKFGZfdNoNPD09DS5PajQx8TEQKVSwdfXt9xs+/fvh4+PDzIyMkrXioqKEBcXZ/bYBmChJyJJGAwqs2+PKzY2FnXr1oWjo6PJfRkZGcjNzQUAdOzYEc7Ozpg4cSJiY2Nx6dIlTJw4EZmZmQgODjb7vbPQE5EUDCXm3x5XRkbGA0946ty5M9atWwcAqFGjBiIjI+Hk5ITg4GAMGDAA2dnZ2LRp0wOP1HkcnNETkRRKBFzrprydqPHx8UY/N2zYECtXrnwqOdjRExEpHDt6IpICr15JRKRw/IYpIiKFq4gTpqoqzuiJiBSOHT0RSYGjGyIihRNxeGVlwUJPRFLgUTdERArHnbFERKRY7OiJSAqc0RMRKRxn9ERECifzjJ6FnoikwNGNIC+dLhL58pXG+bunRUeoNCLcXxQdodIISTsqOkKl8vAv46MHYUdPRFLgjJ6ISOE4uiEiUjiJ98XyhCkiIqVjR09EUuDohohI4bgzlohI4UpEBxCIhZ6IpGCAvB09d8YSESkcO3oikkKJxMdXstATkRRKJB7dsNATkRRkntGz0BORFGQ+6oY7Y4mIFI4dPRFJgaMbIiKFk3l0w0JPRFKQudBzRk9EpHDs6IlICpzRExEpXIm8dZ6FnojkwDNjiYgUTuJL3XBnLBGR0rGjJyIpyHx4JQs9EUmhRMUZPRGRosk8o2ehJyIpyDy64c5YIiKFY0dPRFLgCVNERArHE6aIiBRO5p2xnNETESkcO3oikoLMM3p29ACmLPwAv6X+KjqGMKNGDsa5s0eRnXUFFy/8jNCQEaIjCePm442+26dhRPxXGBrzBV5cNhoOtTWiYwmhtM9FSQXcqirpC32Hrn54MyhQdAxhxoWOxBcRC/HDD4fw+r+GY9u2XViyeCamTwsTHc3iarf0Qr8d01FSWIT/jFmBqEU7UPcfLdF73Qeio1mcEj8Xhgq4VVUqg8EgLL+PxwuiXhoAUM2lGr45ugElJSV4xtNDWJ7zd68JeV2VSoUb12Pwf4ePY1jwuNL19es+Q+9eL+KZZ1tZPFOE+4sWf837+mydCrXGAbv6zYKh+F7/1qB3e7wwewi+778AWdfSLJonJO2oRV/vvsr4uQCAIv1Nsx6/1nOw2RlGJG96rO13796Nf//737hx4wbq1auHkJAQvPzyy2Vuq9Pp8PHHH+PAgQMoKChAly5dEB4ejtq1a5ud+6Ez+lOnTmHnzp3IyspCt27d0L9/f1hbW5fen5WVhXfffRebN282O4ylTft4Aq4n3sC56AsYGRYkOo7FGQwG9H5lILTabKP1goIC2NmpBaUSw666M+p0fA4/TfqqtMgDwNUD0bh6IFpgMsvj56Ji7NmzB9OmTcPkyZPRrVs37N+/Hx988AHc3NzQrl07k+1nzpyJmJgYrFixAmq1GrNmzUJoaCi2bt1qdpZyRzdHjhxBcHAwMjIyUFhYiLlz52Lw4MHQarWl2xQWFiImJsbsIJbW67Ue6PxSJ8wKmw+B/1Mj3IULl5CUdK9TqlmzBoYHD8SQwW/gy1Vfiw1mYbWeqwsrayvkZ2Sh+/IxGB63BiMufYXun42F2sVRdDyLU+LnwpIzeoPBgM8++wyDBw9GUFAQ6tevj7Fjx6JTp044efKkyfapqanYs2cPpk+fjvbt26NVq1ZYtmwZYmJiEB1tfqNRbke/cuVKjB8/HqNHjwYA/P777wgJCUFwcDA2btwIR8eq+QfA1b02pn48EZ/OXIGUZMv+73hl1a1rJ/zff74BAESfOYvln60RnMiyHGrd2+HadfFI3Dh2DgdHLofGyx3+U97Cq16TsOufswEJGwIlfS4suTM1MTERN2/eRJ8+fYzW165dW+b2MTExKCkpgZ+fX+la/fr14eHhgdOnT6N9+/Zm5Sm3o//jjz/wyiuvlP7cpk0bREZGIjk5GePGjUNxcbFZLy7KzGVTce7MBeza8r3oKJXG5StX0b3HvzB0WCiqu2hw6r8/ws3N/NlgVWFle6/nuX3hOo5NXIObv15E3OYjOD5tPdzbNkK9bmLm0qIp6XNhUJl/02q1SE5ONrn975QDAK5duwYA0Ov1eOedd9CxY0e8+eabOHLkSJnZ0tLSUL16dTg4OBitu7m5ISUlxez3Xm6hr1GjBm7eNN4B4u3tjYiICJw6dQrTp0+vcmOPN4cFomXb5lgweQmsra1hbW0N1f+/TvX//rtsbt5Mwc/HT2LLlu/Qp98Q1Kv3LEYMHyQ6lsUU5uYDAJKO/G60fuOn8wCA2s/Xt3imykBJn4uKGN1ERkaiR48eJrfIyEij18rJyQEATJo0Cb169cK6devQuXNnvPvuu/j1V9NDufPz82Fra2uyrlarodfrzX7v5Y5uXnrpJcycObN0bnR/VOPr64v58+dj8uTJSEurWqOPnv26Q1Ndgx/PfGdyX/TNn7F3+37MHD9fQDLLc3HR4NVXX8KJE6dx7dqN0vUrV64iK0sLT886AtNZVtbVVACAtdr4j4SVzb0DD4oKCi2eSRR+Lh4sKCgIgYGmh2NrNMbnWtwv2sHBwfjXv/4FAHjuuedw4cIFrFu3Di+8YHyEn729PQoLTT9jer2+Qkbk5Rb68ePHIyUlBaNHj8aaNWvQuXPn0vv69esHlUqF8PBws0NY0rwPF8PJ2fgX92bQP/HPQX3xdq8RyLz7p6BkYqxd8ylW/3sjwt6fUbrWsUN7uLhocPbsRYHJLCvz8i1ok9Lh3a8Dzn11oHS9foAPACAlKl5UNCGU+LmoiBm9RqMxKepl8fDwAAA0adLEaL1x48Y4fPhwmdtnZWVBp9PBzs6udD09Pb30ucxRbqF3cnLCihUrcPv2bZPZEQD07dsXfn5+OHTokNlBLOX6H0kmaxlpdwAAsWcvWTqOUFlZWny6bBUmTngXOTm5+OmnE2jatBGmTR2PmN/OI3LDDtERLerk/K0I+DIUvdaEIXbzEVT3fgZ+H76JqwejkXE2UXQ8i1Hq58KSQ+bmzZvDyckJ58+fh7+/f+l6QkIC6tWrZ7L9/cMto6Ki0KVLFwDA9evXkZqaCl9fX7PzPNK1bso7YN/d3R1DhgwxOwiJ8dGMRUhOTsHo0UMRNn4U7tzJxLbtuzFz1mLodDrR8Swqcf9p/Dj8U7QPC0SvNWHQZeVhoiT8AAAPzUlEQVQhdstRRC2qmoXNHEr8XFjyWjf29vYYOXIkVq5cCTc3N7Rp0wb79u3DL7/8gvXr1wMAMjIy4OjoCCcnJ7i7u+PVV1/FzJkzsWDBAjg5OWHWrFnw8/ODj4+P2XmkPjO2shB1ZmxlJPLM2MpG1JmxlZW5Z8Z+Vs/8M2PHJz3embFff/01Nm3ahNTUVDRs2BChoaEICAgAADRt2hQhISEIDQ0FAOTl5WHBggU4ePAgDAYDunTpghkzZqBmzZpm52ahrwRY6P/CQv8XFnpj5hb6ZRVQ6N9/zEJfWfAyxUQkhap89UlzsdATkRSq1hk/FYuFnoikwC8eISIixWJHT0RS4IyeiEjhOKMnIlK4EolLPWf0REQKx46eiKTAGT0RkcLJO7hhoSciSbCjJyJSOJ4wRUREisWOnoikIPPhlSz0RCQFecs8Cz0RSYI7Y4mIFE7m0Q13xhIRKRw7eiKSgrz9PAs9EUmCM3oiIoXjjJ6IiBSLHT0RSUHefp6FnogkwRk9EZHCGSTu6VnoiUgKMnf03BlLRKRw7OiJSAoyH17JQk9EUpC3zLPQE5Ek2NETESkcd8YSEZFisaMnIinwOHpB1ti4inz5ysONv4f7OqQdFR2h0jjp5is6gqLIPLphR09EUpC5o+eMnohI4djRE5EUOLohIlK4EoO8oxsWeiKSgrxlnoWeiCQh85mx3BlLRKRw7OiJSAoyH17JQk9EUuBRN0RECifzjJ6FnoikIPPohjtjiYgUjh09EUmBM3oiIoUz8MxYIiJlk3lnLGf0RERP0dWrV+Hj44NvvvnmgdscP34cTZs2NbmdOHGiQjKwoyciKYiY0RcWFmLixInIy8srd7v4+Hg0btwY69evN1p3cXGpkBws9EQkBRGHV65YsQJOTk4P3S4hIQGNGzeGq+vT+bY5jm6ISAolMJh9exynT5/G9u3bsWjRooduGx8fD29v7yd9aw/Fjp6IpFARR91otVpotVqTdY1GA41GY7TdpEmT8NFHH+GZZ54p9zmLioqQmJiIhIQEBAYGIiMjA02aNEFYWBhatWpldmaAhZ6I6JFFRkYiIiLCZD0kJAShoaGlP8+aNQtt2rRB3759H/qc169fh16vR0FBAcLDw2FlZYUNGzZg8ODB+Pbbb9G4cWOzc7PQE5EUKmJnbFBQEAIDA03W/7eb3717N6Kjo/H9998/0nN6e3vj9OnTcHZ2hpXVvWn64sWL0adPH2zYsAFz5841OzcLPRFJoSJ2xv59RFOWb7/9Fnfu3EG3bt2M1ufMmYOvv/4a+/btK/N5/5eVlRUaNWqEW7dumZ0ZYKEnIklY6oSpJUuWoKCgwGitZ8+eCAkJQZ8+fUy2379/P6ZPn45Dhw6VHnVTVFSEuLg49OjRo0IysdATkRQsdQkEd3f3Mtdr1qyJZ599FgCQkZEBR0dHODk5oWPHjnB2dsbEiRMxefJkWFlZYdWqVcjMzERwcHCFZOLhlUREFta5c2esW7cOAFCjRg1ERkbCyckJwcHBGDBgALKzs7Fp06YH/qXxuNjRE5EURF7rJj4+vtyfGzZsiJUrVz6112ehJyIpyPzFIyz0RCSFEokvU8wZPRGRwrGjJyIpyNvPs9ATkSRk/uIRFnoikgILvUSqdXweTb+Z98D7b326Dbc+3WbBROLwd2Fq1MjBCA0dgQZedZF04xZWrYrEioi1omNZlFI/F/zOWInkXUjEpcCpJut1JgyEU+tGuLvnuIBUYvB3YWxc6EgsWTwTS5auxNGjv6JDh3ZYsngmNJpqmL9gueh4FsPPhfJIV+iLs/OQczrOaM0lwBeazq3wx5jFKPjjpqBklsffxV9UKhUmffgetmzdhWnTFwIA/vN/P6NhQy+EvDdcqkKv1M8FRzcSU9mrUW/uKPx5OBqZP/wqOo5QMv8uDAYDer8yEFptttF6QUEB7OzUglJVDkr5XPCEqXLk5+cjISEBTZo0gYODA2JjY7Fx40akpaXB29sbw4YNK71QT1XkPrIf1B61kDAgXHQU4WT/XVy4cKn032vWrIF/vtYbQwa/gc8+XyMwlXhK+VzIPKMv94SpK1euICAgAP3798fLL7+MEydOYNCgQTh37hyqVauGY8eOITAwEFeuXLFU3gqlUtvAffiruLv3OHTXUkXHEYq/i79069oJ6akX8O/VS3DhYjyWfyZvoVfS58LS3xlbmZRb6BcvXoy2bdti9+7d8PX1xdixY/HKK6/ghx9+wGeffYYff/wRnTt3xscff2ypvBWqRt/OsHWrgdQvd4uOIhx/F3+5fOUquvf4F4YOC0V1Fw1O/fdHuLnVFh1LCH4ulKHcQh8VFYWwsDA0a9YMkydPhk6nw9tvvw2VSgUAsLGxwZgxY3DmzBmLhK1oNV7piPz4JOTHXRMdRTj+Lv5y82YKfj5+Elu2fIc+/YagXr1nMWL4INGxhFDS58JgMJh9q6rKLfR2dnbQ6XQAgNq1ayMwMBD29vZG22RnZ8PZ2fnpJXxKVGobuHRtg8x9J0RHEY6/C8DFRYNBg16Hl1ddo/UrV64iK0sLT886gpKJo7TPBUc3D9CpUyfMnz8f169fBwAsXLgQ3t7epffHxMRg5syZePHFF59uyqfAsXkDWNnbITsqVnQU4fi7uGftmk8RNv4do7WOHdrDxUWDs2cvCkoljtI+F4YK+KeqKrfQT548GdnZ2fj8889N7vvhhx8waNAg1K5dGx9++OFTC/i0ODSrDwAoSLghOIl4/F0AWVlafLpsFcaOCcK8uVMQ8NI/EPLecHy7cy1ifjuPyA07REe0OH4ulKPcwyvd3d2xZ88e3L592+Q+f39/7NixAy1btiyd2VclNrWrAwCKsnIEJxGPv4t7PpqxCMnJKRg9eijCxo/CnTuZ2LZ9N2bOWlw6wpSJ0j4XMl+PXmUQuIch2vOfol6aKqkO6adFR6g0Trr5io5QqbRPNu/Inxbu/mZnuJh2yuznEEH6M2OJSA4yd/Qs9EQkhaq8M9Vc/CpBIiKFY0dPRFLg6IaISOFkHt2w0BORFNjRExEpnMwdPXfGEhEpHDt6IpKCwVAiOoIwLPREJIWqfPVJc7HQE5EUqvL15M3FGT0RkcKxoyciKXB0Q0SkcDKPbljoiUgKPGGKiEjheMIUEREpFjt6IpICZ/RERArHo26IiBRO5o6eM3oiIoVjR09EUuDhlURECifz6IaFnoikwJ2xREQKJ3NHz52xREQKx46eiKTAnbFERAon87VuWOiJSArs6ImIFI47Y4mIqEKlpaXhgw8+gL+/P3x8fPDOO+/g8uXLD9w+MzMTEyZMgJ+fH3x9fTFjxgzk5uZWSBYWeiKSgqEC/nnk1zIYMGrUKKSmpmLt2rXYuXMn7O3tMWzYsAcW73HjxiEpKQnr169HREQETpw4gfDw8Ap57yz0RCQFg8Fg9u1R3b59G97e3pg/fz6ef/55eHt7491338Xt27eRkJBgsn1MTAyioqKwcOFCtGjRAv7+/pg3bx727duHW7dumf3eWeiJSAqWLPSurq5YtmwZGjRoAOBe4V+7di3c3NzQpEkTk+2jo6NRq1YtNGrUqHStXbt2UKlUiI6ONvu9c2csEdEj0mq10Gq1JusajQYajabMx0yZMgW7du2CWq3Gl19+CScnJ5Nt0tPT4eHhYbSmVqtRo0YNpKammp1baKFvn7xb5MtTJVQkOgApVqH+ptnPsWLFCkRERJish4SEIDQ0tMzHjBgxAm+//Ta2bNmC9957D5s3b8bzzz9vtE1+fj7UarXJY9VqNXQ6ndm52dETET2ioKAgBAYGmqw/qJsHgMaNGwMA5s+fj7Nnz2Ljxo1YtGiR0Tb29vbQ6/Umj9Xr9XB0dDQzNQs9EdEjK29E87/S09Nx6tQp9OnTByqVCgBgZWWFRo0aIS0tzWR7Dw8PpKenG63p9XpkZmaajHSeBHfGEhFVsJSUFEycOBFnzpwpXSssLERsbCy8vb1Ntvf19UVGRgYSExNL1+7vhG3fvr3ZeVQGmU8XIyJ6CkpKSjBs2DDcvn0bc+bMgUajwapVq/Dzzz9j165dqFOnDu7evYtq1arB3t4eBoMBgwYNQkFBAWbPno2CggJMmzYN7du3x8cff2x2HhZ6IqKnICsrC0uWLMHRo0eRnZ2N9u3bY9KkSWjatCmSk5PRo0cPLFy4EK+//joA4M6dO5g9ezaOHz8OtVqNXr16Ydq0abC3tzc7Cws9EZHCcUZPRKRwLPRERArHQk9EpHDSFvqSkhJ8/vnn6NKlC1q3bo3hw4fj+vXromMJt3r1agwcOFB0DCFycnKwYMECdO/eHT4+Pnj99ddx+PBh0bGEeNxL7FLlJm2h/+KLL7B161bMmzcP27dvh7W1NUaMGFEhpxtXVZs3b8ayZctExxBm6tSpOHbsGObNm4fdu3ejZ8+eCAkJwX//+1/R0SzqSS6xS5WcQUI6nc7Qpk0bw6ZNm0rXsrOzDa1btzbs2rVLYDIxUlNTDaNHjza0adPG0Lt3b8OAAQNER7K49PR0Q5MmTQxHjx41Wh86dKjhgw8+EBNKkPT0dENYWJghMTGxdC0uLs7QpEkTQ0xMjMBk9KSk7Ojj4uKQl5eHDh06lK45OzujefPmFXJJ0Krm4sWLcHJywt69e9G6dWvRcYRwcHDAmjVrTM5CVKlUyMrKEpRKjMe9xC5VflJe6+b+tSbc3d2N1t3c3JCSkiIiklDdu3dH9+7dRccQytnZGf/4xz+M1n7//XecPHkSH330kaBU4j3KJXap8pOyo8/PzwcAk8uCqtXqMq8gR/L5448/EBISgtatW6N///6i4wgzYsQI7Ny5E3369MF7772HCxcuiI5ET0DKQn//lOK/F/WKuiQoVW2nT5/GoEGD4OrqitWrV8PW1lZ0JGEaN26Mli1bYv78+Xj22WexceNG0ZHoCUhZ6J955hkAMLksaHp6usk4h+Syd+9eBAcHo0WLFti4cSOqV68uOpLFpaen4/vvvzf66rzyLrFLlZ+Uhb5Zs2ZwdnZGVFRU6VpOTg5iY2Ph5+cnMBmJ9P3332PSpEl4+eWXsXr1ajg7O4uOJMTjXmKXKj8pd8aq1WoMHjwYy5YtQ+3ateHp6YmlS5fC3d0dPXv2FB2PBEhNTcWMGTPg7++PDz/8EH/++Wfpfba2tlJ19i1btoS/vz/Cw8ONLrH7559/YtiwYaLj0ROQstADwLhx41BcXIzw8HDk5+ejXbt2+Oqrr8r83kZSvkOHDiE/Px8nT55Ely5djO5r27Yttm7dKiiZ5VlZWWHFihVYsmQJwsLCSi+xu3nzZtStW1d0PHoCvEwxEZHCSTmjJyKSCQs9EZHCsdATESkcCz0RkcKx0BMRKRwLPRGRwrHQExEpHAs9EZHCsdATESnc/wNbEpQ3OWBLMAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_cm = pd.DataFrame(confusion_matrix(y,yhat,K) )\n",
    "sn.heatmap(df_cm, annot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {},
   "outputs": [],
   "source": [
    "#least squares multi-class classifier\n",
    "A = np.random.randn(4,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 4, 3, 4]"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([1, 4, 3, 4])"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "[1, 0, 2, 2, 1]"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([1, 0, 2, 2, 1])"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "row_argmax = lambda u:[np.argmax(u[i,:]) for i in range(np.size(u,0))]\n",
    "col_argmax = lambda u:[np.argmax(u[:,i]) for i in range(np.size(u,1))]\n",
    "row_argmax(A)\n",
    "np.argmax(A,1)\n",
    "col_argmax(A)\n",
    "np.argmax(A,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 439,
   "metadata": {},
   "outputs": [],
   "source": [
    "#can use the following to find n-vector predictions of a data set with N examples \n",
    "#stored as nxN 'x', and theta is an nxK \n",
    "fhat = lambda x,theta: row_argmax(x@theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 440,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = 4\n",
    "ycl = np.random.randint(K,size=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 442,
   "metadata": {},
   "outputs": [],
   "source": [
    "#matrix least squares\n",
    "def one_hot(ycl,K):\n",
    "    N = len(ycl)\n",
    "    Y = np.zeros((N,K))\n",
    "    for j in range(K):\n",
    "        Y[np.where(ycl==j),j] = 1\n",
    "    return Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 457,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ls_multiclass(X,ycl,K):\n",
    "    n,N = np.size(X),np.size(X) #debugging code left in on julia companion\n",
    "    theta = npl.lstsq(X.T, (2*one_hot(ycl,K)-1))[0]\n",
    "    yhat = np.argmax(np.matmul(X.T,theta),1)\n",
    "    return theta,yhat\n",
    "\n",
    "setosa = irisD[\"setosa\"]\n",
    "versicolor = irisD[\"versicolor\"]\n",
    "virginica = irisD[\"virginica\"]\n",
    "I1 = np.random.permutation(50)\n",
    "I2 = np.random.permutation(50)\n",
    "I3 = np.random.permutation(50)\n",
    "\n",
    "#4x120 data matrix: 4 features per flower, 120 flowers\n",
    "#transposed shape shows 4 arrays, 1 per feature, with 120 data pts in each\n",
    "Xtrain = np.vstack([setosa[I1[:40],:],\n",
    "          versicolor[I2[:40],:],\n",
    "           virginica[I3[:40],:]]).T \n",
    "#add constant feature array 1, now shape 5x120\n",
    "Xtrain = np.vstack([np.ones((1,120)),Xtrain]) \n",
    "Ytrain = np.vstack([np.full(40,1),np.full(40,2),np.full(40,3)]).ravel()\n",
    "\n",
    "Xtest = np.vstack([setosa[I1[40:],:],\n",
    "          versicolor[I2[40:],:],\n",
    "           virginica[I3[40:],:]]).T \n",
    "Xtest = np.vstack([np.ones((1,30)),Xtest]) \n",
    "Ytest = np.vstack([np.full(10,1),np.full(10,2),np.full(10,3)]).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 458,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([[-1.00000000e+00, -8.73033404e-01,  2.23738964e+00],\n",
       "        [-3.11414326e-16,  1.51403727e-01, -9.91707115e-02],\n",
       "        [ 9.01482844e-16,  4.90757265e-01, -8.49216652e-01],\n",
       "        [ 5.96967799e-16, -4.43634366e-01,  4.73792829e-01],\n",
       "        [-8.24418191e-16, -1.43158749e-01, -1.00380638e+00]]),\n",
       " array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
       "        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2,\n",
       "        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
       "        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
       "        2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2,\n",
       "        2, 2, 2, 0, 2, 2, 2, 2, 2, 1]))"
      ]
     },
     "execution_count": 458,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.,  0.],\n",
       "       [ 0., 40.,  0.],\n",
       "       [ 0.,  0., 40.]])"
      ]
     },
     "execution_count": 458,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.3333333333333333"
      ]
     },
     "execution_count": 458,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0.],\n",
       "       [0., 0., 0.],\n",
       "       [0., 0., 0.]])"
      ]
     },
     "execution_count": 458,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "0.36666666666666664"
      ]
     },
     "execution_count": 458,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta,yhat = ls_multiclass(Xtrain,Ytrain,3)\n",
    "theta,yhat\n",
    "Ctrain = confusion_matrix(Ytrain,yhat,3)\n",
    "Ctrain\n",
    "error_train = error_rate(Ytrain,yhat)\n",
    "error_train\n",
    "yhat = row_argmax(np.matmul(Xtest.T,theta))\n",
    "Ctest = confusion_matrix(Ytest,yhat,3)\n",
    "error_test = error_rate(Ytest,yhat)\n",
    "Ctest\n",
    "error_test\n",
    "\n",
    "#performance is 2x worse than VLMS companion, need to revisit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
